
void initialise (ind * individuals,int n, int intial_col, int xDim, int yDim){

  int x, y;

  //create the grid

 /* grid =  Calloc(xDim,Grid *); //allocate memory for first dimension X

  for (int xx=0; xx<xDim; xx++)
  {
   grid[xx] = Calloc(yDim,Grid);//allocate memory for second dimension Y
   for (int yy=0; yy<yDim; yy++)//set everything to 0
   {
	grid[xx][yy].N=0;
	grid[xx][yy].newN=0;
   }
  } */


  //Initialise the first species (n/2 individuals)
  for (int i = 0; i < (n/2); i++)
  {
   cur_ind[i].x_coord = (int)Rf_runif(0,initial_col);

   cur_ind[i].y_coord = (int)Rf_runif(0,yDim);

   x = cur_ind[i].x_coord;

   y = cur_ind[i].y_coord;

   cur_ind[i].hybrid_index = 0;

   cur_ind[i].mean_negexI = Rf_runif(0,xDim/10);

   cur_ind[i].mean_negexII = Rf_runif(0,xDim/10);

   grid[x][y].N++;
  }

  //Initialise the second species (n/2 individuals)

  for (int i = (n/2); i < n; i++)

  {

   cur_ind[i].x_coord = (int)Rf_runif(xDim - initial_col,initial_col);

   cur_ind[i].y_coord = (int)Rf_runif(0,yDim);

   x = cur_ind[i].x_coord;

   y = cur_ind[i].y_coord;

   cur_ind[i].hybrid_index = 1;

   cur_ind[i].mean_negexI = Rf_runif(0,xDim/10);

   cur_ind[i].mean_negexII = Rf_runif(0,xDim/10);

   grid[x][y].N++;

  }



  for (int xx = 0; xx < xDim; xx++)

  {

   for (int yy = 0; yy < yDim; yy++)

   {

	//allocate 3rd dim of this matrix

	grid[x][y].ind_index = Calloc(grid[x][y].N,int)

	for (int zz = 0; zz < grid[x][y].N; zz++){grid[x][y].ind_index[zz] = -99;}

   }

  }

  for (int i = 0; i < n; i++)

  {

   x = cur_ind[i].x_coord;

   y = cur_ind[i].y_coord;

   grid[x][y].ind_index[first_empty()]= i;

  }



  popsize = n;

}

double choosiness (double diffs, double C, int a, int k){
  return (exp(a-C*k*diffs)/(1+exp(a-C*k*diffs)));
}
//ii is the index when going through the individuals array
//C = "chosiness"
//xx & yy are the individual coordinate
//This function is for 1 individual, so it must be used inside an individual loop
int mate_selection(ind * individuals, int **gd, int ***in, int ii, int xx, int yy, double C)
{
 int z;
 double *HID;
 double *attract;
 double sum_attract = 0;
 double *P_choose;
 HID = Calloc(gd[xx][yy], double); //hybrid index distances
 attract = Calloc(gd[xx][yy], double);//mate attractiveness
 cumPchoose = Calloc(gd[xx][yy], double); //cumulative choosing probabilities

 for (int j = 0; j < gd[xx][yy]; j++)
 {
  int jj = in[xx][yy][j] ;
  HID[j]=individuals[ii].HI - individuals[jj].HI;
  attract[j] = choosiness(HID[j],C,5,100);
  sum_attract += attract[j];
 }
 cumPchoose[0]= 0.0;
 for (int j = 0; j < gd[xx][yy]-1; j++){
   cumPchoose[j] += attract[j] / sum_attract;
 }
 double rn = runif(0.0,sum_attract);
 for (int j = 0; j < gd[xx][yy]-1; j++){
	if (rn <= cumPchoose[j]) return in[xx][yy][j];
 }
 Free(cumPchoose);
 Free(HID);
 Free(attract);
}

// x = hybrid index
double hybrid_fitness(double maxcost, double x, double alpha, double beta){
  return (1-maxcost*Rf_dbeta(x, alpha, beta)/Rf_dbeta(alpha/(alpha+beta), alpha, beta));
}

void reproduce(ind *individuals, int **gd, int ***in, double maxcost,
			   double alpha, double beta, double mut_rate, int xDim)
{
 int x, y;
 int n_offs;
 int nn = 0; //counter for new_inds
 int MI; //mate index
 double E, fitness;
 ind *new_ind = NULL; //define the memory for new indiviuduals

 //Loop through each inds
 for (int i = 0; i < popsize; i++)
 {
  if (individuals[i].alive == 1)
  {
   x = individuals[i].x_coord;
   y = individuals[i].y_coord;
   //calculate individual fitness
   fitness = hybrid_fitness(maxcost,individuals[i].HI,alpha,beta);
   //hassel commins population growth where b=1 (contest competition)
   //generates expected fecundity
   E = fitness*lambda/(1+(lambda-1)*gd[x][y]/K);
   n_offs = Rf_rpois(E);

   if (n_offs > 0)
   {
	//select a mate
	double meanC = (individuals[i].C1+individuals[i].C2)/2;
	MI = mate_selection(individuals, gd, in, i, x, y, meanC);
	//allocate memory for this set of individuals
	new_ind = Realloc(nn+n_offsp, ind);
	//loop through new offspring and initialise them
	for (int j = 0; j < n_offs; j++)
	{
	 new_ind[nn].x_coord = x;
	 new_ind[nn].y_coord = y;
	 new_ind[nn].alive = 1;
	 new_ind[nn].HI = (individuals[i].HI+individuals[MI].hybrid_index)/2;
	 new_ind[nn].mean_negexI =(individuals[i].mean_negexI+individuals[i].mean_negexII)/2;
	 new_ind[nn].mean_negexII =(individuals[MI].mean_negexI+individuals[MI].mean_negexII)/2 ;
	 new_ind[nn].C1 = (individuals[i].C1+individuals[i].C2)/2;
	 new_ind[nn].C2 = (individuals[MI].C1+individuals[MI].C2)/2;

	 mutations(new_ind, nn, mut_rate, xDim);

	 nn++;//increment new_ind counter
	}
   }
  }
 }
 popsize=nn; // adjust popsize

 &cur_ind = &new_ind; //pass pointer
 new_ind = Realloc(0, ind);//allocate zero memory for this set of individuals
}

void mutations(ind* individuals, int id, double mut_rate, int xDim)
{
 double mut;
 //Mutate the 1st allele for dispersal distance
 mut = Rf_rbinom(1,mut_rate);
 if (mut == 1){
	do{
	   individuals[id].mean_negexI = individuals[id].mean_negexI + Rf_runif(-xDim/10, xDim/10);
	}while (individuals[id].mean_negexI < 0.0);
 }
 //Mutate the 2nd allele for dispersal distance
 mut = Rf_rbinom(1,mut_rate);
 if (mut == 1){
	 do{
		individuals[id].mean_negexII = individuals[id].mean_negexII + Rf_runif(-xDim/10, xDim/10);
	 }while (individuals[id].mean_negexII < 0.0);
 }
 //Mutate the 1st allele for choosiness
// mut = Rf_rbinom(1,mut_rate);
// if (mut == 1){
//	 do{
//		individuals[id].C1 = individuals[id].C1 + Rf_runif(-0.1, 0.1);
//	 }while (individuals[id].C1 < 0.0 || individuals[id].C1 > 1.0 );
// }
   //Mutate the 2nd allele for choosiness
// mut = Rf_rbinom(1,mut_rate);
// if (mut == 1){
//	 do{
//		individuals[id].C2 = individuals[id].C2 + Rf_runif(-0.1, 0.1);
//	 }while (individuals[id].C2 < 0.0 || individuals[id].C2 > 1.0 );
// }
}

//mort = disoersal mortality
void disperse(ind *individuals, int xDim, int yDim, double mort)
{
 double x_rand, y_rand;
 double R1, dist, meanDist, angle;
 int new_x, new_y;
 int alive;

 int n = popsize;
 //at the moment there is no dispersal probability - every individual moves.
 //Loop through each individual
 for (int i = 0; i < n; i++){
	alive = Rf_rbinom(1,mort);
	if (alive == 1)
	{
	 x = individuals[i].x_coord;
	 y = individuals[i].y_coord;
	 meanDist = (individuals[i].mean_negexI + individuals[i].mean_negexII)/2;
	 //Calculate destination
	 x_rand = Rf_runif(0.0,1.0);
	 y_rand = Rf_runif(0.0,1.0);
	 //Reflective boundaries on the x
	 do{
	  R1 = Rf_runif(0.0000001,1.0);
	  dist = (-1.0*meanDist)*log(R1);
	  angle = Rf_runif(0.0,1.0)*2.0*pi;
	  new_x=(int)(dist*cos(angle)+x_rand+x);
	  new_y=(int)(dist*sin(angle)+y_rand+y);
	 } while (new_x < 0 || new_x >=(xDim-1));
	 //Wrapped boundaries on the y
	 if (new_y >= yDim) {
	   do{
		 new_y = new_y - yDim;
	   }while(new_y < yDim);
	 }
	 individuals[i].x_coord = new_x;
	 individuals[i].y_coord = new_y;
	}
	else
	{
	 individuals[i].alive = 0;
	 popsize--;
    }
 }
}

